c----------------------------------------------------------------------------
c Program: 2D_Potential_FMM - A fast multipole boundary element
c Method (BEM) code for analyzing large-scale, general 2D
c potential problems (governed by Laplace equation)
c using constant elements.
c
c Developers: Dr. Naoshi Nishimura at Kyoto University, Japan;
c Dr. Yijun Liu at the University of Cincinnati, Cincinnati, OH.
c
c Version: V.1.20.
c Released: October 1, 2008.
c
c Copyright(c)2004-2008 By Kyoto University and University of Cincinnati.
c This code is intended for educational use only. No part of
c the code can be used for any commercial applications/
c distributions without prior written permissions of the
c original developers.
c
c----------------------------------------------------------------------------

      implicit real*8(a-h,o-z)
      integer, allocatable:: ia(:)
      complex*16, allocatable :: am(:)
      character*80 Prob_Title
      call CPU_Time(time0)
      open (4, file='input.fmm',    status='old')
      open (5, file='input.dat',    status='old')
      open (3, file='output.dat',   status='unknown')
      open (7, file='phi_boundary.plt', status='unknown')
      open (8, file='xy.plt',       status='unknown')
      open (9, file='phi_domain.plt',     status='unknown')

c Input the parameters
      read(4,*)   maxl, levmx, nexp, ntylr, tolerance
      read(4,*)   maxia, ncellmx, nleafmx, mxl, nwksz
      read(5,'(a80)')   Prob_Title
      read(5,*)   n, nfield
      write(3,'(a80)') Prob_Title
      write(*,'(a80)') Prob_Title

c Estimate the maximum numbers of the cells and leaves,
c and size of the preconditioning matrix, etc.
      if(ncellmx.le.0) ncellmx = max(4*n/maxl,100)
      if(nleafmx.le.0) nleafmx = max(ncellmx/2,100)
      if(nwksz.le.0)
      nwksz = maxl*maxl*nleafmx
      ligw = 20
      lrgw = 1+n*(mxl+6)+mxl*(mxl+3)
      iwksz = n+3*nleafmx+1
      allocate (ia(maxia))

c Load the addresses (pointers) associated with the locations of the
c variables to be stored in the large array "am"
      call lpointer(lp, ln, maxia, ia, n, nexp, ntylr, ncellmx,
     & levmx, ligw, lrgw, nwksz, iwksz, nfield,
     & l_n, l_x, l_y, l_node, l_dnorm,
     & l_bc, l_a, l_b, l_xmax,
     & l_xmin, l_ymax, l_ymin, l_ielem, l_itree,
     & l_level, l_loct, l_numt, l_ifath, l_lowlev,
     & l_maxl, l_levmx, l_nexp, l_ntylr, l_tolerance,
     & l_ncellmx, l_nleafmx, l_mxl, l_u, l_ax,
     & l_sb, l_sx, l_ligw, l_lrgw, l_igwk,
     & l_rgwk, l_nwksz, l_iwksz, l_rwork, l_iwork,
     & l_xfield, l_nfield, l_f)

c Estimate the memory usage

      maxa = lp
      write(3,100) maxa*16/1.D6
      write(*,100) maxa*16/1.D6
  100 format(/'Memory size of the large block am =', f12.1,' Mb'/)

c Allocate the large block 'am'
      allocate (am(maxa))

c Assign the parameters to the array am()
      call assigni(n, am(l_n ))
      call assigni(maxl, am(l_maxl ))
      call assigni(levmx, am(l_levmx ))
      call assigni(nexp, am(l_nexp ))
      call assigni(ntylr, am(l_ntylr ))
      call assignd(tolerance, am(l_tolerance))
      call assigni(ncellmx, am(l_ncellmx))
      call assigni(nleafmx, am(l_nleafmx))
      call assigni(mxl, am(l_mxl ))
      call assigni(ligw, am(l_ligw ))
      call assigni(lrgw, am(l_lrgw ))
      call assigni(nwksz, am(l_nwksz ))
      call assigni(iwksz, am(l_iwksz ))
      call assigni(nfield, am(l_nfield))

c Call the FMM BEM main program
      call fmmmain(maxa, maxia, am, ia,
     & am(l_n), am(l_x), am(l_y), am(l_node),
     & am(l_dnorm), am(l_bc), am(l_a), am(l_b),
     & am(l_xmax), am(l_xmin), am(l_ymax), am(l_ymin),
     & am(l_ielem), am(l_itree), am(l_level), am(l_loct),
     & am(l_numt), am(l_ifath), am(l_lowlev), am(l_maxl),
     & am(l_levmx), am(l_nexp), am(l_ntylr), am(l_tolerance),
     & am(l_ncellmx), am(l_nleafmx), am(l_mxl), am(l_u),
     & am(l_ax), am(l_nfield), am(l_xfield), am(l_f),
     & am(l_sb), am(l_sx), am(l_igwk), am(l_rgwk),
     & am(l_ligw), am(l_lrgw), am(l_nwksz), am(l_iwksz),
     & am(l_rwork), am(l_iwork))

c Estimate the total CPU time
      call CPU_Time(time)
      write(3,*)
      write(*,*)
      write(3,*) ' Total CPU time used =', time-time0, '(sec)'
      write(*,*) ' Total CPU time used =', time-time0, '(sec)'
      stop
      end

c----------------------------------------------------------------------------
c Definition of Variables:
c
c maxa = maximum size of the array am
c maxia = maximum number of variables allowed
c am = a large array storing the variables for the SLATEC GMRES solver
c ia = an array storing the locations of the variables in the array am
c
c n = number of elements (= number of nodes)
c x = coordinates of the nodes
c y = coordinates of the end points of the elements
c node = element connectivity
c dnorm = normal at each node
c bc = BC type and value
c
c a = multipole expansion moments
c b = local expansion coefficients
c xmax,xmin = maximum and minimum x coordinate
c ymax,ymin = maximum and minimum y coordinate
c ielem = ielem(i) gives the original element number for i-th element in
c the quad-tree structure c itree = itree(c) gives the cell location of c-th cell within each
c tree level
c loct = elements included in the c-th cell are listed starting at
c the loct(c)-th place in the array ielem
c numt = numt(c) gives the number of elements included in the c-th cell
c ifath = ifath(c) gives the number of the parent cell of the c-th cell
c level = level l cells start at the level(l)-th cell in the tree
c lowlev = number of the tree levels
c
c maxl = maximum number of elements allowed in a leaf
c levmx = maximum number of levels allowed in the tree structure
c nexp = number of terms in multipole expansion
c ntylr = number of terms in local expansion
c tolerance = GMRES solution convergence tolerance
c ncellmx = maximum number of cells allowed in the tree
c nleafmx = maximum number of leaves allowed in the tree
c mxl = maximum dimension of Krylov subspace (used in GMRES)
c
c u = first stores b vector; then solution vector of system Ax = b
c ax = resulting vector of multiplication Ax
c nfield = number of the field points inside the domain
c xfield = coordinates of the field points inside the domain
c f = values of the potential at the field points inside the domain
c
c The following variables and arrays are used in the SLATEC GMRES solver:
c sb,sx,igwk,rgwk,ligw,lrgw,nwksz,iwksz,rwork,iwork
c
c----------------------------------------------------------------------------

      subroutine lpointer(lp, ln, maxia, ia, n, nexp, ntylr, ncellmx,
     & levmx, ligw, lrgw, nwksz, iwksz, nfield,
     & l_n, l_x, l_y, l_node, l_dnorm,
     & l_bc, l_a, l_b, l_xmax,
     & l_xmin, l_ymax, l_ymin, l_ielem, l_itree,
     & l_level, l_loct, l_numt, l_ifath, l_lowlev,
     & l_maxl, l_levmx, l_nexp, l_ntylr, l_tolerance,
     & l_ncellmx, l_nleafmx, l_mxl, l_u, l_ax,
     & l_sb, l_sx, l_ligw, l_lrgw, l_igwk,
     & l_rgwk, l_nwksz, l_iwksz, l_rwork, l_iwork,
     & l_xfield, l_nfield, l_f)

         dimension ia(maxia)
         lp= 1

         l_n = l_address( 1,maxia,ia,lp,4,1)
         l_x = l_address( 2,maxia,ia,lp,8,n*2)
         l_y = l_address( 3,maxia,ia,lp,8,n*2)
         l_node = l_address( 4,maxia,ia,lp,4,n*2)
         l_dnorm = l_address( 5,maxia,ia,lp,8,n*2)
         l_bc = l_address( 6,maxia,ia,lp,8,n*2)
         l_a = l_address( 7,maxia,ia,lp,16,(nexp+1)*ncellmx)
         l_b = l_address( 8,maxia,ia,lp,16,(ntylr+1)*ncellmx)
         l_xmax = l_address( 9,maxia,ia,lp,8,1)
         l_xmin = l_address(10,maxia,ia,lp,8,1)
         l_ymax = l_address(11,maxia,ia,lp,8,1)
         l_ymin = l_address(12,maxia,ia,lp,8,1)
         l_ielem = l_address(13,maxia,ia,lp,4,n)
         l_itree = l_address(14,maxia,ia,lp,4,ncellmx)
         l_level = l_address(15,maxia,ia,lp,4,levmx+1)
         l_loct = l_address(16,maxia,ia,lp,4,ncellmx)
         l_numt = l_address(17,maxia,ia,lp,4,ncellmx)
         l_ifath = l_address(18,maxia,ia,lp,4,ncellmx)
         l_lowlev = l_address(19,maxia,ia,lp,4,1)
         l_maxl = l_address(20,maxia,ia,lp,4,1)
         l_levmx = l_address(21,maxia,ia,lp,4,1)
         l_nexp = l_address(22,maxia,ia,lp,4,1)
         l_ntylr = l_address(23,maxia,ia,lp,4,1)
         l_tolerance = l_address(24,maxia,ia,lp,8,1)
         l_ncellmx = l_address(25,maxia,ia,lp,4,1)
         l_nleafmx = l_address(26,maxia,ia,lp,4,1)
         l_mxl = l_address(27,maxia,ia,lp,4,1)
         l_u = l_address(28,maxia,ia,lp,8,n)
         l_ax = l_address(29,maxia,ia,lp,8,n)
         l_sb = l_address(30,maxia,ia,lp,8,n)
         l_sx = l_address(31,maxia,ia,lp,8,n)
         l_ligw = l_address(32,maxia,ia,lp,4,1)
         l_lrgw = l_address(33,maxia,ia,lp,4,1)
         l_igwk = l_address(34,maxia,ia,lp,4,ligw)
         l_rgwk = l_address(35,maxia,ia,lp,8,lrgw)

         l_nwksz = l_address(36,maxia,ia,lp,4,1)
         l_iwksz = l_address(37,maxia,ia,lp,4,1)
         l_rwork = l_address(38,maxia,ia,lp,8,nwksz)
         l_iwork = l_address(39,maxia,ia,lp,4,iwksz)
         l_xfield = l_address(40,maxia,ia,lp,8,nfield*2)
         l_nfield = l_address(41,maxia,ia,lp,4,1)
         l_f = l_address(42,maxia,ia,lp,8,nfield)
         return
      end

c----------------------------------------------------------------------------

      integer function l_address(ln,maxia,ia,lp,ibyte,length)
         dimension ia(maxia)
         l_address = lp
         ia(ln) = lp

         iu = 16
         inc = (ibyte*length-1)/iu+1
         lp = lp+inc
         if(ln .gt. maxia) then
            write(*,*)'!Specified # of variables maxia',maxia,'is too small'
            stop
         endif
         return
      end

c----------------------------------------------------------------------------
      subroutine assigni(i,ii)
         integer i,ii
         ii = i
         return
      end

      subroutine assignd(d,dd)
         real*8 d,dd
         dd = d
         return
      end

c----------------------------------------------------------------------------
      subroutine fmmmain(maxa, maxia, am, ia, n,x,y,node,dnorm,bc,
     & a,b, xmax,xmin,ymax,ymin,ielem,itree,level,loct,numt,
     & ifath,lowlev,maxl,levmx,nexp,ntylr,tolerance,ncellmx,
     & nleafmx,mxl,u,ax,nfield,xfield,f,sb,sx,igwk,rgwk,
     & ligw,lrgw,nwksz,iwksz,rwork,iwork)

         implicit real*8(a-h,o-z)
         complex*16 am(maxa), a,b
         dimension ia(maxia),ja(1),a(0:nexp,ncellmx),b(0:ntylr,ncellmx),
     &    x(2,n),y(2,n),node(2,n),dnorm(2,n),bc(2,n),
     &    ielem(n),itree(ncellmx),level(0:levmx),loct(ncellmx),
     &    numt(ncellmx),ifath(ncellmx), u(n),ax(n),sb(n),sx(n),
     &    igwk(ligw),rgwk(lrgw),rwork(nwksz),iwork(iwksz),
     &    xfield(2,nfield),f(nfield)

         external matvec, msolve

c Input parameters and prepare the BEM model
         call prep_model(n,x,y,node,bc,dnorm,xfield,nfield,maxl,levmx,
     &    nexp,ntylr,tolerance,xmin,xmax,ymin,ymax)

c Generate the quad-tree structure for the elements
         call tree(n,x,xmax,xmin,ymax,ymin,ielem,itree,level,loct,numt,
     &    ifath,lowlev,maxl,levmx,ncellmx,nleafmx,nwksz,iwork)

c Compute the right-hand-side vector b with the FMM
         call fmmbvector(n,x,y,node,dnorm,bc,u,ax,a,b,xmax,xmin,ymax,ymin,
     &    ielem,itree,level,loct,numt,ifath,
     &    nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)

c Solve the BEM system of equations Ax=b with the fast multipole BEM

c Prepare parameters for calling the iterative solver GMRES
c (SLATEC GMRES solver is used, which is available at www.netlib.org.
c See the documentation for the SLATEC GMRES solver for more information
c about the following related parameters)

         nelt = 1
         isym = 0
         itol = 0
         tol = tolerance
         iunit = 3
         igwk(1) = mxl
         igwk(2) = mxl
         igwk(3) = 0
         igwk(4) = 1
         igwk(5) = 10
         do i=1,n
            ax(i) = 0.d0
         enddo

         write(*,*) 'Call Equation Solver GMRES ...'
         call dgmres(n,u,ax, nelt,ia,ja,am,isym, matvec,msolve,itol,tol,
     &    itmax,iter,er,ierr,iunit,sb,sx,rgwk,lrgw,igwk,ligw,
     &    rwork,iwork)
         write(3,*) ' Error indicator from GMRES:', ierr
         write(*,*) ' Error indicator from GMRES:', ierr

c Output the boundary solution
         do i=1,n
            u(ielem(i)) = ax(i)
         enddo

         write(3,*) ' Fast Multipole BEM Solution:'
         do i=1,n
            write(3,*) i, u(i)
            write(7,*) i, u(i)
         enddo

c Evaluate the field inside the domain and output the results
         call domain_field(nfield,xfield,f,n,x,y,bc,node,dnorm,u)
         return
      end
c----------------------------------------------------------------------------

      subroutine prep_model(n,x,y,node,bc,dnorm,xfield,nfield,maxl,
     & levmx,nexp,ntylr,tolerance,
     & xmin,xmax,ymin,ymax)

         implicit real*8(a-h,o-z)
         dimension x(2,*),y(2,*),node(2,*),bc(2,*),dnorm(2,*),xfield(2,*)
         write(*,2) n, maxl, levmx, nexp, tolerance
         write(3,2) n, maxl, levmx, nexp, tolerance
    2    format(' Total number of elements =', I12
         & /' Max. number of elements in a leaf =', I12
         & /' Max. number of tree levels =', I12
         & /' Number of terms used in expansions =', I12
         & /' Tolerance for convergence =', D12.3)
         write(*,*)
         write(3,*)

c Input the mesh data
         read(5,*)
         do i=1,n
            read(5,*) itemp, y(1,i), y(2,i)
         enddo
         read(5,*)
         do i=1,n
            read(5,*) itemp, node(1,i), node(2,i), bc(1,i), bc(2,i)
         enddo

c Input the field points inside the domain
         if (nfield .gt. 0) then
            read(5,*)
            do i=1,nfield
               read(5,*) itemp, xfield(1,i), xfield(2,i)
            enddo
         endif

c Compute mid-nodes and normals of the elements
         do i=1,n
            x(1,i) = (y(1,node(1,i)) + y(1,node(2,i)))*0.5
            x(2,i) = (y(2,node(1,i)) + y(2,node(2,i)))*0.5
            h1 = y(2,node(2,i)) - y(2,node(1,i))
            h2 = -y(1,node(2,i)) + y(1,node(1,i))
            el = sqrt(h1**2 + h2**2)
            dnorm(1,i) = h1/el
            dnorm(2,i) = h2/el
         enddo

c Determine the square bounding the problem domain (Largest cell used in FMM)
         xmin=x(1,1)
         xmax=x(1,1)
         ymin=x(2,1)
         ymax=x(2,1)

         do 10 i=2,n
            if(x(1,i).le.xmin) then
               xmin=x(1,i)
            elseif(x(1,i).ge.xmax) then
               xmax=x(1,i)
            endif
            if(x(2,i).le.ymin) then
               ymin=x(2,i)
            elseif(x(2,i).ge.ymax) then
               ymax=x(2,i)
            endif
   10    continue

         scale = 1.05d0 ! Make the square slightly larger
         xyd = max(xmax-xmin,ymax-ymin)/2.d0
         xyd = xyd*scale
         cx = (xmin+xmax)/2.d0
         cy = (ymin+ymax)/2.d0
         xmin = cx-xyd
         xmax = cx+xyd
         ymin = cy-xyd
         ymax = cy+xyd

c Output nodal coordinates for plotting
         do i = 1,n
            write(8,*) x(1,i), x(2,i)
         enddo
         return
      end

c----------------------------------------------------------------------------
      subroutine tree(n,x,xmax,xmin,ymax,ymin,ielem,itree,level,loct,
         & numt,ifath,lowlev,maxl,levmx,ncellmx,nleafmx,
         & nwksz,iwork)

         implicit real*8(a-h,o-z)
         complex*16 a,b
         dimension x(2,*),ielem(*),itree(*),level(0:*),loct(*),numt(*),
         & ifath(*), iwork(*), nwk(4)

         do i=1,n
            ielem(i) = i ! Store the original element numbers in ielem
         enddo

c For the level 0 cell (largest cell)
         itree(1) = 0
         level(0) = 1
         level(1) = 2
         loct(1) = 1
         ifath(1) = 1
         numt(1) = n
         ndivx = 1
         lowlev = 1
         nleaf = 0
         nswa = 0

c For cells on level 1 to the lowest level (leaves)
         do 10 lev=1,levmx
            levp = lev-1
            levn = lev+1
            level(levn) = level(lev)
            if(level(lev).eq.level(levp)) goto 900
            ndivxp = ndivx
            ndivx = 2*ndivxp

            dxp = (xmax-xmin)/ndivxp ! Parent cell size
            dyp = (ymax-ymin)/ndivxp do 11 inp=level(levp),level(lev)-1
            itrp = itree(inp)
            if(numt(inp).gt.maxl .or.
     &       (lev.le.2 .and. numt(inp).ne.0) ) then
               itrpx = mod(itrp,ndivxp)
               itrpy = itrp/ndivxp
               xsep = xmin+(itrpx + 0.5d0)*dxp
               ysep = ymin+(itrpy + 0.5d0)*dyp

               call bisec(x,ielem(loct(inp)),numt(inp),ysep,nsepy, 2)
               call bisec(x,ielem(loct(inp)),nsepy-1,
               xsep,nsepx1,1)
               call bisec(x,ielem(loct(inp)+nsepy-1),
     &          numt(inp)-nsepy+1,xsep,nsepx2,1)

               nwk(1) = nsepx1-1
               nwk(2) = nsepy-nsepx1
               nwk(3) = nsepx2-1
               nwk(4) = numt(inp)-nsepy-nsepx2+2
               locc = loct(inp)

               do 12 icldy=0,1
                  do 12 icldx=0,1
                     icld = icldy*2+icldx+1

                     if(nwk(icld).gt.0) then
                        nrel = level(levn)
                        if(nrel.gt.ncellmx) then
                           write(*,*) " ncellmx error"
                           stop
                        endif
                        itree(nrel) = ((itrpy*2+icldy)*ndivxp + itrpx)*2
                        & +icldx

                        loct(nrel) = locc
                        numt(nrel) = nwk(icld)
                        ifath(nrel) = inp
                        lowlev=lev

c Leaves:

                        if((lev.ne.1) .and.
     &                   (numt(nrel).le.maxl .or. lev.eq.levmx)) then
                           nleaf = nleaf+1
                           if(nleaf.gt.nleafmx) then
                              write(*,*) " nleafmx error"
                              stop
                           endif
                           nleaf3 = nleaf*3 - 1
                           iwork(nleaf3) = nrel ! Store cell number (icell)
                           iwork(nleaf3+1) = nswa + 1 ! Location of pre-cond'er
                           iwork(nleaf3+2) = 1 ! Initial value of switch isw
                           nswa = nswa + numt(nrel)**2
                           if(nswa.gt.nwksz) then
                              write(*,*) " nwksz error"
                              stop
                           endif
                        endif

                        level(levn) = nrel + 1
                        locc = locc + nwk(icld)
                     endif
   12          continue
            endif
   11       continue
   10    continue

  900    iwork(1) = nleaf ! Store number of leaves in iwork(1)
         write(3,*) ' Number of tree levels =', lowlev
         write(*,*) ' Number of tree levels =', lowlev
         write(3,*) ' Number of leaves =', nleaf
         write(*,*) ' Number of leaves =', nleaf
         write(3,*) ' Number of cells =', nrel
         write(*,*) ' Number of cells =', nrel
         write(3,*)
         write(*,*)
         return
      end

c----------------------------------------------------------------------------
      subroutine bisec(x,ielem,n,xsep,nsep,ic)
         implicit real*8(a-h,o-z)
         dimension x(2,*),ielem(*)
         nsep = 1
         if(n.le.0) return

         do ifr=1,n
            if(x(ic,ielem(ifr)).le.xsep) then
               if(ifr.ne.nsep) then
                  istore = ielem(nsep)
                  ielem(nsep) = ielem(ifr)
                  ielem(ifr) = istore
               endif
               nsep = nsep + 1
            endif
         enddo

         return
      end

c----------------------------------------------------------------------------
      subroutine fmmbvector(n,x,y,node,dnorm,bc,u,ax,a,b,xmax,xmin,
     & ymax,ymin,ielem,itree,level,loct,numt,ifath,
     & nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)
         implicit real*8(a-h,o-z)
         complex*16 a,b
         dimension a(0:nexp,ncellmx),b(0:ntylr,ncellmx),
     &    x(2,*),y(2,*),node(2,*),dnorm(2,*),bc(2,*),u(*),ax(*),
     &    ielem(*),itree(*),level(0:*),loct(*),numt(*),ifath(*),
     &    rwork(*),iwork(*)

c Switch the BC type
         do i=1,n
            if(bc(1,i) .eq. 1.) then
               bc(1,i) = 2.d0
            else
               bc(1,i) = 1.d0
            endif
         enddo

         do i=1,n
            u(i)= bc(2,ielem(i))
            ax(i) = 0.d0
         enddo

c Apply the FMM to conpute the right-hand side vector b

         call upward(u,n,y,node,dnorm,bc,a,xmax,xmin,ymax,ymin,ielem,
         & itree,level,loct,numt,ifath,nexp,ncellmx,lowlev,maxl)

         call dwnwrd(u,ax,n,x,y,node,dnorm,bc,a,b,xmax,xmin,ymax,ymin,
         & ielem,itree,level,loct,numt,ifath,nexp,ntylr,ncellmx,
         & lowlev,maxl,rwork,iwork)

c Store b vector in u and switch the BC type back
         do i=1,n
            u(i) = - ax(i)
            if(bc(1,i) .eq. 1.) then
               bc(1,i) = 2.d0
            else
               bc(1,i) = 1.d0
            endif
         enddo
         return
      end

c----------------------------------------------------------------------------
      subroutine matvec(n,u,ax,nelt,ia,ja,am,isym)
         implicit
         real*8(a-h,o-z)
         complex*16 am, a, b
         dimension u(*),ax(*),ia(*),ja(*),am(*)

c Retrieve the pointers
         l_n = ia(1)
         l_x = ia(2)
         l_y = ia(3)
         l_node = ia(4)
         l_dnorm = ia(5)
         l_bc = ia(6)
         l_a = ia(7)
         l_b = ia(8)
         l_xmax = ia(9)
         l_xmin = ia(10)
         l_ymax = ia(11)
         l_ymin = ia(12)
         l_ielem = ia(13)
         l_itree = ia(14)
         l_level = ia(15)
         l_loct = ia(16)
         l_numt = ia(17)
         l_ifath = ia(18)
         l_lowlev = ia(19)
         l_maxl = ia(20)
         l_levmx = ia(21)
         l_nexp = ia(22)
         l_ntylr = ia(23)
         l_tolerance = ia(24)
         l_ncellmx = ia(25)
         l_nleafmx = ia(26)
         l_mxl = ia(27)
         l_u = ia(28)
         l_ax = ia(29)
         l_sb = ia(30)
         l_sx = ia(31)
         l_ligw = ia(32)
         l_lrgw = ia(33)
         l_igwk = ia(34)
         l_rgwk = ia(35)
         l_nwksz = ia(36)
         l_iwksz = ia(37)
         l_rwork = ia(38)
         l_iwork = ia(39)

c Evaluate matrix-vector multiplication Ax using the fast multipole BEM

         call upward(u, am(l_n), am(l_y), am(l_node),
     &    am(l_dnorm), am(l_bc), am(l_a), am(l_xmax),
     &    am(l_xmin), am(l_ymax), am(l_ymin), am(l_ielem),
     &    am(l_itree), am(l_level), am(l_loct), am(l_numt),
     &    am(l_ifath), am(l_nexp), am(l_ncellmx), am(l_lowlev),
     &    am(l_maxl))

         call dwnwrd(u,ax, am(l_n), am(l_x), am(l_y),
     &    am(l_node), am(l_dnorm), am(l_bc), am(l_a),
     &    am(l_b), am(l_xmax), am(l_xmin), am(l_ymax),
     &    am(l_ymin), am(l_ielem), am(l_itree), am(l_level),
     &    am(l_loct), am(l_numt), am(l_ifath), am(l_nexp),
     &    am(l_ntylr), am(l_ncellmx), am(l_lowlev), am(l_maxl),
     &    am(l_rwork), am(l_iwork))

         return
      end

c----------------------------------------------------------------------------
      subroutine msolve(n,r,z,nelt,ia,ja,am,isym,rwork,iwork)
         implicit real*8(a-h,o-z)
         complex*16 am(*)
         dimension r(*),z(*),ia(*),ja(1),rwork(*),iwork(*)

c Load the pointers
         l_loct = ia(16)
         l_numt = ia(17)

c Compute the preconditioning matrix
         call msolveinv(r,z,rwork,iwork, am(l_loct), am(l_numt))
         return
      end

c----------------------------------------------------------------------------
      subroutine msolveinv(r,z,rwork,iwork,loct,numt)
         implicit real*8(a-h,o-z)
         dimension r(*),z(*),iwork(*),rwork(*),loct(*),numt(*)

         nleaf = iwork(1)
         do l = 1,nleaf
            l3 = l*3-1
            inod = iwork(l3)

            indr = iwork(l3+1)
            indx = loct(inod)
            indi = indx+3*nleaf+1
            nr = numt(inod)

            call dcopy(nr,r(indx),1,z(indx),1)
            call dluax(rwork(indr),nr,nr,z(indx),iwork(l3+2),
     &       iwork(indi),icon)

            if(icon.ne.0) then
               write(*,*) " dluax error, icon =", icon
               stop
            endif
            iwork(l3+2) = 2
         enddo

         return
      end

c----------------------------------------------------------------------------
c This subroutine solves linear system of equations Ax=b by LU decomposition.
c
c a .... given regular coefficient matrix.
c k .... given adjustable dimension for array a.
c n .... given order of matrix a.
c b .... given constant vector.
c isw .... given control information: c if 1, solve equations entirely.
c if 2, solve equations with last LU-decomposed entries.
c ip .... auxiliary 1 dimensioned array, size is n.
c transposition vector which represents
c row-exchanging by partial pivoting.
c icon.... resultant condition code.
c
c Slave subroutines used (available at www.netlib.org):
c dgetrf, dgetrs
c
c----------------------------------------------------------------------------

      subroutine dluax(a,k,n,b,isw,ip,icon)
         implicit real*8 (a-h,o-z)
         dimension a(k,n),b(n),ip(n)
         data ione/1/

         icon = 30000
         if(isw.eq.1) go to 1000
         if(isw.eq.2) go to 1100
         go to 8000
 1000    call dgetrf(n,n,a,k,ip,icon)
         isw = 2
 1100    call dgetrs('n',n,ione,a,k,ip,b,n,icon)
 8000    continue

         return
      end

c----------------------------------------------------------------------------
      subroutine upward(u,n,y,node,dnorm,bc,a,xmax,xmin,ymax,ymin,ielem,
     & itree,level,loct,numt,ifath,nexp,ncellmx,lowlev,maxl)

         implicit real*8(a-h,o-z)
         complex*16 a,b, z0,zi
         dimension a(0:nexp,ncellmx),
     &    y(2,*),node(2,*),dnorm(2,*),bc(2,*),u(*),
     &    ielem(*),itree(*),level(0:*),loct(*),numt(*),ifath(*)

         do i=1,level(lowlev+1)-1
            do k=0,nexp
               a(k,i) = (0.d0,0.d0) ! Clear multipole moments
            enddo
         enddo

         do 10 lev=lowlev,2,-1 ! Loop from leaf to level 2 cells (Upward)
            ndivx = 2**lev
            dx = (xmax-xmin)/ndivx ! Determine cell size
            dy = (ymax-ymin)/ndivx

            do 20 icell=level(lev),level(lev+1)-1 ! Loop for level l cells
               itr = itree(icell)
               itrx = mod(itr,ndivx)
               itry = itr/ndivx ! Position of the cell
               cx = xmin+(itrx + 0.5d0)*dx
               cy = ymin+(itry + 0.5d0)*dy ! Center of the cell

c Multipole expansion
               if(numt(icell).le.maxl .or. lev.eq.lowlev) then ! Compute moment
                  call moment(a(0,icell),y,node,ielem(loct(icell)),
                  & numt(icell),nexp,cx,cy,u(loct(icell)),
                  & bc,dnorm)
               endif

c M2M translation
               if(lev.ne.2) then ! Do M2M translation to form moments
                  cxp = xmin+(int(itrx/2)*2 + 1)*dx
                  cyp = ymin+(int(itry/2)*2 + 1)*dy ! Center of parent cell
                  z0 = cmplx(cx-cxp, cy-cyp) ! (z_c - z_c')
                  io = ifath(icell) ! Cell no. of parent cell
                  zi = (1.d0,0.d0)
                  do k=0,nexp
                     do m=k,nexp
                        a(m,io) = a(m,io) + zi*a(m-k,icell) ! Use M2M
                     enddo
                     zi = zi*z0/(k+1)
                  enddo
               endif

   20       continue
   10    continue
         return
      end

c----------------------------------------------------------------------------
      subroutine moment(a,y,node,ielem,num,nexp,cx,cy,u,bc,dnorm)
         implicit real*8(a-h,o--y)
         implicit complex*16(z)
         complex*16 a(0:*)
         dimension y(2,*),node(2,*),ielem(*),u(*), bc(2,*),dnorm(2,*)

         do i=1,num ! Over elements in the leaf
            nelm = ielem(i) ! Element number
            n1 = node(1,nelm) ! Two ends of the element
            n2 = node(2,nelm)
            z1 = cmplx(y(1,n1)-cx, y(2,n1)-cy)
            z2 = cmplx(y(1,n2)-cx, y(2,n2)-cy)
            zwbar = conjg(z2 - z1)
            zwbar = zwbar/abs(zwbar) ! omega bar
            zp1 = z1*zwbar
            zp2 = z2*zwbar
            znorm = cmplx(dnorm(1,nelm),dnorm(2,nelm)) ! complex normal n
            if(bc(1,nelm) .eq. 1.d0) then ! Assign values to phi and q
               phi = 0.D0
               q = u(i)
            else if(bc(1,nelm) .eq. 2.d0) then
               phi = u(i)
               q = 0.D0
            endif

c Compute moments:
            a(0) = a(0) - (zp2-zp1)*q ! G kernel
            do k=1,nexp
               a(k) = a(k) + (zp2-zp1)*znorm*phi ! F kernel
               zp1 = zp1*z1/(k+1)
               zp2 = zp2*z2/(k+1)
               a(k) = a(k) - (zp2-zp1)*q ! G kernel
            enddo
         enddo
         return
      end

c----------------------------------------------------------------------------
      subroutine dwnwrd(u,ax,n,x,y,node,dnorm,bc,a,b,xmax,xmin,
     & ymax,ymin,ielem,itree,level,loct,numt,ifath,
     & nexp,ntylr,ncellmx,lowlev,maxl,rwork,iwork)

         implicit real*8(a-h,o-z)
         complex*16 a,b, z0,zi,zo,zp

         dimension a(0:nexp,ncellmx),b(0:ntylr,ncellmx),x(2,*),y(2,*),
     &    node(2,*),dnorm(2,*),bc(2,*),u(*),ax(*),ielem(*),
     &    itree(*),level(0:*),loct(*),numt(*),ifath(*),
     &    rwork(*),iwork(*)

         data pi/3.141592653589793D0/
         pi2 = pi*2.d0
         do i=1,level(lowlev+1)-1
            do k=0,ntylr
               b(k,i) = (0.d0,0.d0)
            enddo
         enddo
         do i=1,n
            ax(i) = 0.d0
         enddo

         leaf = 0
         indr = 1
         indi = 1
         do 110 lev=2,lowlev ! Downward from level 2 cells to leaf cells
            ndivx = 2**lev
            dx = (xmax-xmin)/ndivx
            dy = (ymax-ymin)/ndivx
            do 120 icell=level(lev),level(lev+1)-1 ! Loop for level l cells
               itr = itree(icell)
               itrx = mod(itr,ndivx)
               itry = itr/ndivx ! Position of the cell
               cx = xmin+(itrx + 0.5d0)*dx
               cy = ymin+(itry + 0.5d0)*dy ! Center of the cell
               itrxp = itrx/2
               itryp = itry/2

c From the parent cell (use L2L)
               if(lev.ne.2) then
                  cxp = xmin+(itrxp*2+1)*dx
                  cyp = ymin+(itryp*2+1)*dy ! Center of the parent cell
                  z0 = cmplx(cx-cxp, cy-cyp)
                  io = ifath(icell) ! Cell no. of the parent cell
                  zi = (1.d0,0.d0)
                  do k=0,ntylr
                     do m=0,ntylr-k

                        b(m,icell) = b(m,icell) + zi*b(k+m,io) ! L2L translation
                     enddo
                     zi = zi*z0/(k+1)
                  enddo
               endif

               do 130 jcell=level(lev),level(lev+1)-1
                  jtr = itree(jcell)
                  jtrx = mod(jtr,ndivx)
                  jtry = jtr/ndivx
                  jtrxp = jtrx/2
                  jtryp = jtry/2

c The parents must be neighbours
                  if(iabs(itrxp-jtrxp).gt.1 .or. iabs(itryp-jtryp).gt.1)
     &             goto 130

c For non-neighbours (cells in interaction list) (use M2L)
                  if(iabs(itrx-jtrx).gt.1 .or. iabs(itry-jtry).gt.1) then
                     ccx = xmin + (jtrx + 0.5d0)*dx
                     ccy = ymin + (jtry + 0.5d0)*dy ! Center of the j cell
                     z0 = cmplx(cx-ccx, cy-ccy)
                     b(0,icell) = b(0,icell) - log(z0)*a(0,jcell)
                     zo = 1.
                     do m=1,nexp+ntylr
                        zo = zo/z0
                        kmin = max(0,m-nexp)
                        kmax = min(m,ntylr)
                        sgn = (-1.0)**kmin
                        do k=kmin,kmax
                           b(k,icell) = b(k,icell) + sgn*zo*a(m-k,jcell) ! M2L
                           sgn = -sgn
                        enddo
                        zo = zo*m
                     enddo

c Contribution from neighbouring leaves (use direct)
                     elseif(numt(jcell).le.maxl .or.
                     & numt(icell).le.maxl .or. lev.eq.lowlev) then
                     if(icell.eq.jcell) then
                        leaf = leaf+1
                        leaf3 = leaf*3-1
                        if(iwork(leaf3).ne.icell) then
                           write(3,*)
                           leaf,iwork(1),iwork(leaf3),icell,'check'
                           write(*,*) " icell error"
                           stop
                        endif
                        indr = iwork(leaf3+1)
                        indi = iwork(leaf3+2)
                     endif

                     call direct(ielem(loct(icell)),ielem(loct(jcell)),
     &                node,x,y,numt(icell),numt(jcell),dnorm,
     &                ax(loct(icell)),u(loct(jcell)),icell,jcell,
     &                rwork(indr),indi,bc) ! Direct integration
                  endif
  130          continue

c Compute Ax if reach a leaf (Evaluate local expansion at each
c collocation point)
               if(numt(icell).le.maxl .or. lev.eq.lowlev) then
                  fact = 1.d0
                  do itylr=1,ntylr
                     fact = fact/itylr
                     b(itylr,icell) = b(itylr,icell)*fact
                  enddo

                  do in=1,numt(icell)
                     inax = loct(icell) + in-1 ! Element number in the tree
                     indx = ielem(inax) ! Original element number
                     zp = b(ntylr,icell)
                     z0 = cmplx(x(1,indx)-cx, x(2,indx)-cy)
                     do itylr=ntylr-1,0,-1
                        zp = zp*z0 + b(itylr,icell) ! Local expansion
                     enddo
                     zp = zp/pi2
                     ax(inax) = ax(inax) + dreal(zp) ! Array Ax
                  enddo
               endif

  120       continue
  110    continue

         return
      end

c----------------------------------------------------------------------------
      subroutine direct(inod,jnod,node,x,y,ni,nj,dnorm,ax,u,icell,jcell,
     & amat,isw,bc)

     implicit real*8(a-h,o-z)
         dimension inod(*),jnod(*),node(2,*),x(2,*),y(2,*),
     &    dnorm(2,*),ax(*),u(*),amat(ni,*), bc(2,*)
         data pi/3.141592653589793D0/
         pi2 = pi*2.d0
         do j = 1, nj

            jind = jnod(j)
            al = sqrt((y(1,node(1,jind))- y(1,node(2,jind)))**2 +
            (y(2,node(1,jind)) - y(2,node(2,jind)))**2) ! Element length

            do i = 1, ni
               iind = inod(i)
               x11 = y(1,node(1,jind)) - x(1,iind)
               x21 = y(2,node(1,jind)) - x(2,iind)
               x12 = y(1,node(2,jind)) - x(1,iind)
               x22 = y(2,node(2,jind)) - x(2,iind)
               r1 = sqrt(x11**2 + x21**2)
               r2 = sqrt(x12**2 + x22**2)
               d = x11*dnorm(1,jind) + x21*dnorm(2,jind)
               t1 = -x11*dnorm(2,jind) + x21*dnorm(1,jind)
               t2 = -x12*dnorm(2,jind) + x22*dnorm(1,jind)
               ds = abs(d)
               dtheta = datan2(ds*al,ds**2+t1*t2)
               aa = (-dtheta*ds + al + t1*log(r1) - t2*log(r2))/pi2
               if(d .lt. 0.d0) dtheta = -dtheta

               bb = -dtheta/pi2
               if(iind .eq. jind) bb = 0.5d0
               if(bc(1,jind) .eq. 1.) then
                  ax(i) = ax(i) - aa*u(j)
                  if(icell.eq.jcell .and. isw.eq.1) then ! Store coefficients in
                     amat(i,j) = - aa ! first iteration
                  endif
               else if(bc(1,jind) .eq. 2.) then
                  ax(i) = ax(i) + bb*u(j)
                  if(icell.eq.jcell .and. isw.eq.1) then ! Store coefficients in
                     amat(i,j) = bb ! first iteration
                  endif
               endif
            enddo
         enddo
         return
      end

c----------------------------------------------------------------------------
      subroutine domain_field(nfield,xfield,f,n,x,y,bc,node,dnorm,u)
         implicit real*8(a-h,o-z)

         dimension xfield(2,*), f(*), x(2,*),y(2,*),bc(2,*),node(2,*),
         & dnorm(2,*),u(*)
         data pi/3.141592653589793D0/
         pi2 = pi*2.d0
         do i=1,nfield
            f(i) = 0.d0
         enddo
         do j=1,n ! Loop over all elements
            if(bc(1,j).eq.1) then
               f0 = bc(2,j)
               df0 = u(j)

            else if(bc(1,j).eq.2) then

               f0 = u(j)
               df0 = bc(2,j)
            endif
            al= sqrt((y(1,node(2,j)) - y(1,node(1,j)))**2 +
            (y(2,node(2,j)) - y(2,node(1,j)))**2) ! Element length

            do i=1,nfield ! Loop over all field points inside the domain
               x11 = y(1,node(1,j)) - xfield(1,i)
               x21 = y(2,node(1,j)) - xfield(2,i)
               x12 = y(1,node(2,j)) - xfield(1,i)
               x22 = y(2,node(2,j)) - xfield(2,i)
               r1 = sqrt(x11**2 + x21**2)
               r2 = sqrt(x12**2 + x22**2)
               d = x11*dnorm(1,j) + x21*dnorm(2,j)
               t1 = -x11*dnorm(2,j) + x21*dnorm(1,j)
               t2 = -x12*dnorm(2,j) + x22*dnorm(1,j)
               ds = abs(d)
               theta1 = datan2(t1,ds)
               theta2 = datan2(t2,ds)
               dtheta = theta2 - theta1

               aa = (-dtheta*ds + al + t1*log(r1) - t2*log(r2))/pi2
               if(d .lt. 0.d0) dtheta = -dtheta
               bb = -dtheta/pi2
               f(i) = f(i) + aa*df0 - bb*f0
            enddo
         enddo

c Output results
         do i=1,nfield
            write(9,20) xfield(1,i), f(i)
         enddo
   20    format(1x, 4E18.8)
         return
      end

c----------------------------------------------------------------------------
